\documentclass[a4paper,twoside]{ctexart}
\usepackage{geometry}
\geometry{margin=1cm,vmargin={0pt,1cm}}
\setlength{\topmargin}{-2cm}
\setlength{\paperheight}{23cm}
\setlength{\paperwidth}{18cm}
\setlength{\textheight}{19.6cm}
\setlength{\textwidth}{15cm}
\usepackage{makecell}
\usepackage{fancyhdr}
\usepackage{siunitx}
\usepackage{amssymb}
\usepackage{indentfirst}
\setlength{\parindent}{0.5em}

\pagenumbering{arabic}

% useful packages.
\usepackage{multirow}
\usepackage{caption}
\usepackage{mathrsfs}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{enumerate}
\usepackage{xcolor,graphicx,float,subfigure}
\usepackage{epstopdf}
\usepackage{multicol}
\usepackage{fancyhdr}
\usepackage{layout}
\usepackage{listings}
\usepackage{diagbox}
\lstset{language=Matlab}
\lstset{breaklines}
\lstset{extendedchars=false}
\usepackage[colorlinks,linkcolor=blue]{hyperref}
\usepackage{xcolor}
\usepackage{cite}
\usepackage[numbers,sort&compress]{natbib}
\setcitestyle{open={},close={}}
%\usepackage{natbibspacing}
%\renewcommand{\refname}{}
\usepackage{anyfontsize}

% English theorem environment
\theoremstyle{definition}
\newtheorem{proposition}[section]{Proposition}
\newtheorem{corollary}[section]{Corollary}
\newtheorem{definition}{Definition}[section]
\newtheorem{remark}[section]{Remark}
\newtheorem{example}[definition]{Example}
\newtheorem{exercise}[definition]{Exercise}
\newtheorem{theorem}[definition]{Theorem}
\newtheorem{lemma}[definition]{Lemma}
\newenvironment{solution}{\begin{proof}[Solution]}{\end{proof}}
\renewcommand{\proofname}{\noindent Proof}
% some common command
\newcommand{\dif}{\mathrm{d}}
\newcommand{\avg}[1]{\left\langle #1 \right\rangle}
\newcommand{\pdfrac}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\op}{\odot}
\newcommand{\Eabs}{E_{\mathrm{abs}}}
\newcommand{\Erel}{E_{\mathrm{rel}}}
\newcommand{\Ediv}{\mathrm{div}}%\div是除号
\newcommand{\lrq}[1]{\left( #1 \right)}
\newcommand{\avint}[1]{\frac{1}{\left|#1\right|}\int_{#1}}

\newcommand{\upcite}[1]{\textsuperscript{\textsuperscript{\cite{#1}}}} 

\makeatletter
\newcommand\sixteen{\@setfontsize\sixteen{17pt}{6}}
\renewcommand{\maketitle}{\bgroup\setlength{\parindent}{0pt}
\begin{flushleft}
\sixteen\bfseries \@title
\medskip
\end{flushleft}
\texttt{\@author}
\egroup}
\renewcommand{\maketag@@@}[1]{\hbox{\m@th\normalsize\normalfont#1}}
\makeatother
\setlength{\abovecaptionskip}{0.cm}
\setlength{\floatsep}{0.cm}
\setlength{\intextsep}{0.cm}

%\CTEXsetup[format={\Large\bfseries}]{section}

\title{微分方程数值解\ 第6次作业}
\author{李之琪 22235056}

\begin{document}
\maketitle
\section*{coding parts}
\subsection*{part 1}
对于对流扩散方程，我们在空间上分别实现二阶、四阶算法，并比较达到相同精度所
需要的点数和用时。由于扩散项的存在，若想使用$k = O(h)$的方法，只能
使用隐式方法。这里我们选用Crank-Nicholson方法。

取$a = -1, \eta = 1$，初值函数$f(x) = \sin{x}$，终止时刻$T = 1$。采用加细网格的方法计算近似的
精确解。下面给出为了达到给定精度所需的点数和消耗的CPU时间。
\newline
\begin{table}[!htp]
  \centering
  \begin{tabular}{c|cc}
    \hline
    Order&2&4\\
    \hline
     $M_P$&11&12\\
    \hline
    CPU time(s)&1.88e$-$2&2.74e$-$2\\
    \hline
  \end{tabular}
  \caption{$\varepsilon = 1e-2$}
  \label{table:1}
\end{table}
\newline
\begin{table}[!htp]
  \centering
  \begin{tabular}{c|cc}
    \hline
    Order&2&4\\
    \hline
     $M_P$&35&37\\
    \hline
    CPU time(s)&2.94e$-$1&3.06e$-$1\\
    \hline
  \end{tabular}
  \caption{$\varepsilon = 1e-3$}
  \label{table:2}
\end{table}\\
这里CPU时间的测量采用连续十次测量取平均的方式得到。我们观察到，相较
于二阶方法，四阶方法并不能有效地减少所需的点数，但需要付出更高的代价。
(主要计算量来自线性方程组的求解，因此随空间阶数提高计算时间的增长有
限)因此，我们认为对于对流扩散方程，二阶方法就足够了。

以上结果可以通过运行
\texttt{CrankNicholson\_diffusionConvection\_res1.m}和\\
\texttt{CrankNicholson\_diffusionConvection\_res2.m}分别得到。
\subsection*{part 2}
BDF法的RAS如下图所示，注意RAS是图中的无界区域。

\begin{figure}[!htp]   
  \centering
  \includegraphics[width=12cm]{Pictures/F0.eps}
  \caption{RAS of BDF, $3 \le p \le 6$}
  \label{fig:1}
\end{figure}

该图可以通过运行\texttt{plot\_BDF\_RAS.m}得到。
\section*{2.2.1}
对Simpson方法$\dfrac{u_{j}^{n+1}-u_j^{n-1}}{2k} =
\dfrac{1}{6}(Q_4u_j^{n+1}+4Q_4u_j^{n}+Q_4u_j^{n-1})$，代入$u_j^n =
\dfrac{1}{\sqrt{2\pi}}\hat{v}^n(\omega)e^{i\omega x}$。注意到
\begin{eqnarray*}
  \begin{aligned}
    D_0u_j^{n} &=
    \frac{1}{\sqrt{2\pi}}\hat{v}^n(\omega)\frac{i\sin{\xi}}{h}e^{i\omega x},\\
    D_0D_+D_-u_j^{n} &=
    \frac{-4}{\sqrt{2\pi}}\hat{v}^n(\omega)\frac{i\sin{\xi }\sin^2{\frac{\xi}{2}}}{h^3}e^{i\omega x},
  \end{aligned}
\end{eqnarray*}
这里$\xi = \omega h$。结合$Q_4 = D_0(I - \dfrac{h^2}{6}D_+D_-)$，可得
\begin{eqnarray*}
  \begin{aligned}
    \hat{v}^{n+1}(\omega) - \hat{v}^{n-1}(\omega) &=
    \frac{\sigma }{3}i\left[\frac{3\sin{\xi} + 2\sin{\xi}\sin^2{\frac{\xi}{2}}}{3}\hat{v}^{n+1}(\omega)
    + \frac{4(3\sin{\xi} + 2\sin{\xi}\sin^2{\frac{\xi}{2}})}{3}\hat{v}^{n}(\omega)\right.\\
  &+ \left.\frac{3\sin{\xi} + 2\sin{\xi}\sin^2{\frac{\xi}{2}}}{3}\hat{v}^{n-1}(\omega)\right],
  \end{aligned}
\end{eqnarray*}
这里$\sigma = \dfrac{k}{h}$。整理得
\begin{eqnarray*}
  \begin{aligned}
    \left[1-\frac{\sigma}{9}\left(3\sin{\xi} +
        2\sin{\xi}\sin^2{\frac{\xi}{2}}\right)i\right]\hat{v}^{n+1}(\omega)&-\frac{4\sigma}{9}\left(3\sin{\xi}
      +
      2\sin{\xi}\sin^2{\frac{\xi}{2}}\right)i\hat{v}^{n}(\omega)\\
    &-\left[1+\frac{\sigma}{9}\left(3\sin{\xi}
        +
        2\sin{\xi}\sin^2{\frac{\xi}{2}}\right)i\right]\hat{v}^{n-1}(\omega)
    = 0,
  \end{aligned}
\end{eqnarray*}
对应的特征方程为
\begin{eqnarray*}
  \begin{aligned}
    \left[1-\frac{\sigma}{9}\left(3\sin{\xi} +
        2\sin{\xi}\sin^2{\frac{\xi}{2}}\right)i\right]z^2&-\frac{4\sigma}{9}\left(3\sin{\xi}
      + 2\sin{\xi}\sin^2{\frac{\xi}{2}}\right)iz\\
    &-\left[1+\frac{\sigma}{9}\left(3\sin{\xi}
        + 2\sin{\xi}\sin^2{\frac{\xi}{2}}\right)i\right] = 0,
  \end{aligned}
\end{eqnarray*}
解得
\begin{eqnarray*}
  \begin{aligned}
    z_{1,2} = \frac{2\sigma t i \pm \sqrt{1-3\sigma^2 t^2}}{1-\sigma t i},
  \end{aligned}
\end{eqnarray*}
这里$t(\xi) = \dfrac{3\sin{\xi} + 2\sin{\xi}\sin^2{\frac{\xi}{2}}}{9}
\in (-\dfrac{1}{2},\dfrac{1}{2})$。我们断言：若存在某个$t(\xi)$使得
$1-3\sigma^2 t^2 < 0$，则必有$|z_1| > 1$或$|z_2| > 1$。这是因为韦达定
理告诉我们,
\begin{eqnarray*}
  \begin{aligned}
    |z_1z_2| &= \left|\frac{1+\sigma t i}{1-\sigma t i}\right| = 1,\\
    |z_1+z_2| &= \left|\frac{4\sigma t i}{1-\sigma t i}\right| = \left|\frac{4\sigma t
    i - 4\sigma^2 t^2}{1 +
      \sigma^2 t^2}\right|.
  \end{aligned}
\end{eqnarray*}
第一式说明，$|z_{1,2}| \le 1$当且仅当$|z_{1,2}| = 1$。但第二式给出
\begin{eqnarray*}
  \begin{aligned}
    |z_1+z_2|^2 &= \frac{16\sigma^2 t^2 + 16\sigma^4 t^4}{(1 +
      \sigma^2 t^2)^2} > 4.
  \end{aligned}
\end{eqnarray*}
最后一个不等式等价于$(3\sigma^2 t^2-1)(\sigma^2t^2+1) > 0$。这与
$|z_{1,2}| = 1$矛盾，从而说明了上面的断言。


若对任意$t(\xi)$，$1-3\sigma^2 t^2 \ge 0$，令$\sigma^2 t^2 = s \le
\dfrac{1}{3}$，注意到
\begin{eqnarray*}
  \begin{aligned}
    |z_{1,2}|^2 = \left|\frac{2 \sigma t i \pm \sqrt{1-3\sigma^2
          t^2}}{1-\sigma t i}\right|^2 =
   \frac{[s(2\pm \sqrt{1-3s})]^2 + (-2s \pm \sqrt{1-3s})^2}{(1+s)^2}
   \le 1.
  \end{aligned}
\end{eqnarray*}
结合$t(\xi) \in (-\dfrac{1}{2},\dfrac{1}{2})$，我们得到当$\sigma <
\dfrac{2\sqrt{3}}{3}$时，该方法稳定。
\end{document}

